Celem poniższej analizy było znalezienie przyczyny malającej długości śledzi wyławianych w Europie. Analiza przeprowadzona została na podstawie zbioru danych z pomiarami śledzi i warunków w jakich żyli w ciągu ostatnich 60 lat. Po wyczyszczeniu danych i przeprowadzeniu analizy stwierdzono, że największy wpływ na rozmiar śledzi miała temperatura przy powierzchni wody. Im wyższa temperatura tym mniejszy śledź.
W raporcie wykorzystano następujące biblioteki:
library(knitr)
library(DT)
library(dplyr)
library(tidyr)
library(ggplot2)
library(plotly)
library(ggcorrplot)
library(caret)
Dane zawierają pomiary dotyczące wielkości śledzi i warunków w jakich żyją z ostatnich 60 lat.
Dane zostały zebrane z połowów komercyjnych jednostek.
df <- read.csv("./data/sledzie.csv", na.strings="?")
df <- tbl_df(df)
| Nazwa atrybutu | Opis |
|---|---|
| X | numer obserwacji |
| length | długość złowionego śledzia [cm] |
| cfin1 | dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1] |
| cfin2 | dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2] |
| chel1 | dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1] |
| chel2 | dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2] |
| lcop1 | dostępność planktonu [zagęszczenie widłonogów gat. 1] |
| lcop2 | dostępność planktonu [zagęszczenie widłonogów gat. 2] |
| fbar | natężenie połowów w regionie [ułamek pozostawionego narybku] |
| recr | roczny narybek [liczba śledzi] |
| cumf | łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku] |
| totaln | łączna liczba ryb złowionych w ramach połowu [liczba śledzi] |
| sst | temperatura przy powierzchni wody [°C] |
| sal | poziom zasolenia wody [Knudsen ppt] |
| xmonth | miesiąc połowu [numer miesiąca] |
| nao | oscylacja północnoatlantycka [mb] |
Zbiór danych zawiera 52582 pomiarów. W zbiorze znajduje się 10094 niepełnych obserwacji, co stanowi 19.2 % wszystkich pomiarów. Liczba rekordów z wartościami NA jest zbyt duża, aby je usunąć.
Liczba brakujących wartości w poszczególnych kolumnach przedstawia się następująco:
kable(df %>% summarise_all(funs(sum(is.na(.)))))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 1581 | 1536 | 1555 | 1556 | 1653 | 1591 | 0 | 0 | 0 | 0 | 1584 | 0 | 0 | 0 |
Ze względu na to, iż wartości atrybutów często występują w podobnych grupach, to znaczy, że wartości atrybutów w sąsiadujących obserwacjach często są takie same. Brakujące wartości uzupełniono na podstawie sąsiednich rekordów - górnego lub dolnego.
na_columns <- colnames(df)[colSums(is.na(df)) > 0]
cleaned_df <- df %>% fill(na_columns, .direction="updown")
Wyczyszczony zbiór danych składa się z 52582 wierszy (obserwacji) i 16 kolumn (atrybutów).
kable(summary(cleaned_df %>% select(X:lcop2)))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | |
|---|---|---|---|---|---|---|---|---|
| Min. : 0 | Min. :19.0 | Min. : 0.0000 | Min. : 0.0000 | Min. : 0.000 | Min. : 5.238 | Min. : 0.3074 | Min. : 7.849 | |
| 1st Qu.:13145 | 1st Qu.:24.0 | 1st Qu.: 0.0000 | 1st Qu.: 0.2778 | 1st Qu.: 2.469 | 1st Qu.:13.427 | 1st Qu.: 2.5479 | 1st Qu.:17.808 | |
| Median :26291 | Median :25.5 | Median : 0.1111 | Median : 0.7012 | Median : 5.750 | Median :21.435 | Median : 7.0000 | Median :24.859 | |
| Mean :26291 | Mean :25.3 | Mean : 0.4457 | Mean : 2.0255 | Mean :10.003 | Mean :21.215 | Mean : 12.8079 | Mean :28.419 | |
| 3rd Qu.:39436 | 3rd Qu.:26.5 | 3rd Qu.: 0.3333 | 3rd Qu.: 1.7936 | 3rd Qu.:11.500 | 3rd Qu.:27.193 | 3rd Qu.: 21.2315 | 3rd Qu.:37.232 | |
| Max. :52581 | Max. :32.5 | Max. :37.6667 | Max. :19.3958 | Max. :75.000 | Max. :57.706 | Max. :115.5833 | Max. :68.736 |
kable(summary(cleaned_df %>% select(fbar:nao)))
| fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|
| Min. :0.0680 | Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. :12.77 | Min. :35.40 | Min. : 1.000 | Min. :-4.89000 | |
| 1st Qu.:0.2270 | 1st Qu.: 360061 | 1st Qu.:0.14809 | 1st Qu.: 306068 | 1st Qu.:13.60 | 1st Qu.:35.51 | 1st Qu.: 5.000 | 1st Qu.:-1.89000 | |
| Median :0.3320 | Median : 421391 | Median :0.23191 | Median : 539558 | Median :13.86 | Median :35.51 | Median : 8.000 | Median : 0.20000 | |
| Mean :0.3304 | Mean : 520367 | Mean :0.22981 | Mean : 514973 | Mean :13.88 | Mean :35.51 | Mean : 7.258 | Mean :-0.09236 | |
| 3rd Qu.:0.4560 | 3rd Qu.: 724151 | 3rd Qu.:0.29803 | 3rd Qu.: 730351 | 3rd Qu.:14.16 | 3rd Qu.:35.52 | 3rd Qu.: 9.000 | 3rd Qu.: 1.63000 | |
| Max. :0.8490 | Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :14.73 | Max. :35.61 | Max. :12.000 | Max. : 5.08000 |
p <- ggplot(gather(cleaned_df %>% select(-X, -xmonth)), aes(x = value)) +
geom_histogram() +
facet_wrap(~key, scales="free") +
theme_light()
ggplotly(p)
Poniżej zaprezentowano wykresy zmiany poszczególnych atrybutów na przestrzeni kolejnych obserwacji. Zmianę atrybutów w czasie zaprezentowano za pomocą linii trendu.
plankton_df <- cleaned_df %>%
select(X, cfin1:lcop2) %>%
gather(plankton_name, plankton_density, cfin1:lcop2)
p <- ggplot(plankton_df, aes(x=X, y=plankton_density, color = plankton_name)) +
geom_smooth() +
labs(title = "Zmiana dostępności planktonu", x="Nr obserwacji",
y="Zagęszczenie planktonu") +
scale_color_discrete(name=NULL) +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=fbar)) +
geom_smooth() +
labs(title = "Zmiana natężenia połowów w regionie", x="Nr obserwacji",
y="Natężenie połowów w regionie [ułamek pozostawionego narybku]") +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=recr)) +
geom_smooth() +
labs(title = "Zmiana rocznego narybku", x="Nr obserwacji",
y="Roczny narybek [liczba śledzi]") +
scale_y_continuous(labels=scales::comma) +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=cumf)) +
geom_smooth() +
labs(title = "Zmiana łącznego rocznego natężenia połowów w regionie",
x="Nr obserwacji", y="Ułamek pozostawionego narybku") +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=totaln)) +
geom_smooth() +
labs(title = "Zmiana łącznej liczby ryb złowionych w ramach połowu",
x="Nr obserwacji", y="Liczba śledzi") +
scale_y_continuous(labels=scales::comma) +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=sst)) +
geom_smooth() +
labs(title = "Zmiana temperatury przy powierzchni wody", x="Nr obserwacji",
y="Temperatura [°C]") +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=sal)) +
geom_smooth() +
labs(title = "Zmiana poziomu zasolenia wody", x="Nr obserwacji",
y="Zasolenie [Knudsen ppt]") +
theme_light()
ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=nao)) +
geom_smooth() +
labs(title = "Zmiana oscylacji północnoatlantyckiej", x="Nr obserwacji",
y="Oscylacja północnoatlantycka [mb]") +
theme_light()
ggplotly(p)
Zmianę wielkości śledzia zaprezentowano na interaktywnym wykresie za pomocą linii trendu. Oś x oznacza liczbę porządkową - numer obserwacji, natomiast oś y rozmiar śledzia w cm.
p <- ggplot(cleaned_df, aes(x=X, y=length)) +
geom_smooth() +
labs(x="Nr obserwacji", y="Długość [cm]") +
theme_light()
ggplotly(p)
Korelacje pomiędzy poszczególnymi atrybutami przedstawiono na interaktywnej mapie korelacji. Do obliczenia współczynnika korelacji wykorzystano metodę Pearsona.
Aby wyświetlić nazwy skorelowanych zmiennych i wartości korelacji należy najechać kursorem na poszczególne komórki na wykresie.
cor_matrix <- cleaned_df %>%
select(-X) %>%
cor(use = "all.obs", method="pearson")
corr_plot <- cor_matrix %>%
ggcorrplot(type="lower", legend.title = "Współczynnik Pearsona") +
labs(x = 'Atrybut 1', y = 'Atrybut 2') +
theme_light() +
theme(axis.ticks = element_blank())
ggplotly(corr_plot)
Na podstawie wykresu prezentującego korelację pomiędzy atrybutami można zauważyć, że atrybuty lcop1 i chel1 są ze sobą najsilniej skorelowane. Atrybut length jest najsiliniej dodatnio skorelowany z atrybutami fbar oraz lcop1, natomiast najsilniej skorelowany ujemnie z atrybutami sst oraz nao.
W celu redukcji korelacji między atrybutami wykorzystano funkcję ‘findCorrelation’ z pakietu ‘caret’ z parametrem cutoff równym 0.8, która zwraca atrybuty do usunięcia.
attributes_to_remove <- cor_matrix %>% findCorrelation(cutoff = 0.8, names = TRUE)
Atrybuty wybrane do usunięcia: lcop2, cumf, chel1.
Do predykcji pominięto atrybuty lcop2, cumf, chel1 oraz X. Zbiór danych podzielono na zbiór uczący i testowy w proporcjach 75/25. Zbiór walidacyjny został utworzony przy użyciu wielokrotnej oceny krzyżowej (ang. repeated cross-validation) z liczbą podziałów równą 2 i liczbą powtórzeń równą 5.
in_training_data <- createDataPartition(y = cleaned_df$length, p = 0.75, list = FALSE)
training_data <- cleaned_df[in_training_data, ] %>% select(-c(X, attributes_to_remove))
testing_data <- cleaned_df[-in_training_data, ]
ctrl <- trainControl(method = "repeatedcv", number = 2, repeats = 5)
Poniższy wykres przedstawia podobieństwo rozkładów danych treningowych i testowych.
ggplot() +
geom_density(aes(length, fill = "Treningowy"), training_data, alpha = 0.6) +
geom_density(aes(length, fill = "Testowy"), testing_data, alpha = 0.6) +
labs(x = "Długość [cm]", y = "Gęstość", fill = "Zbiór danych") +
theme_light()
Do przewidywania długości śledzia wykorzystano metodę Random Forest. Do oceny dokładności predykcji wykorzystano miary R2 i RMSE.
fit <- train(length ~ .,
data = training_data,
method = "rf",
trControl = ctrl,
ntree = 10)
fit
## Random Forest
##
## 39438 samples
## 11 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 19720, 19718, 19719, 19719, 19719, 19719, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 1.166928 0.5021753 0.9225294
## 6 1.148083 0.5181004 0.9045890
## 11 1.151060 0.5160672 0.9054773
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 6.
Wartości miar dla zbioru testowego:
prediction <- predict(fit, testing_data)
post_resample <- postResample(pred = prediction,
obs = testing_data$length)
post_resample
## RMSE Rsquared MAE
## 1.1372567 0.5251405 0.8957069
Poniższy wykres przedstawia wartości zbioru testowego oraz wartości przewidziane przez regresor.
prediction_comparison_df <- tibble(X = testing_data$X,
actual = testing_data$length,
predicted = prediction)
ggplot(prediction_comparison_df, aes(x = X)) +
geom_smooth(aes(y = actual, color = "Wartość rzeczywista")) +
geom_smooth(aes(y = predicted, color = "Wartość przewidziana")) +
labs(color = "Wartości", x = "Nr obserwacji", y = "Długość [cm]") +
theme_light()
ggplot(varImp(fit)) +
labs(x = "Atrybut", y = "Ważność") +
theme_light()
Patrząc na powyższy wykres okazało się, że w przewidywaniu rozmiaru śledzia najważniejsza była temperatura przy powierzchni wody (sst). Mniejszy wpływ miały również natężenie połowów w regionie (fbar) oraz roczny narybek (recr).
Na podstawie całej analizy można stwierdzić, że na rozmiar śledzia największy wpływ miała temperatura przy powierzchni wody. Im wyższa temperatura tym mniejszy śledź.